Sector-coupling with pypsa#
Note
If you have not yet set up Python on your computer, you can execute this tutorial in your browser via Google Colab. Click on the rocket in the top right corner and launch “Colab”. If that doesn’t work download the .ipynb file and import it in Google Colab.
Then install the following packages by executing the following command in a Jupyter cell at the top of the notebook.
!pip install pypsa pandas numpy matplotlib plotly highspy
import pypsa
import pandas as pd
import numpy as np
import plotly.io as pio
import plotly.offline as py
import matplotlib.pyplot as plt
pd.options.plotting.backend = "plotly"
Previous electricity-only PyPSA model#
To explore sector-coupling options with PyPSA, let’s load the capacity expansion model we built for the electricity system and add sector-coupling technologies and demands on top.
Some of the sector-coupling technologies we are going to add have multiple ouputs (e.g. CHP plants producing heat and power). PyPSA can automatically handle links have more than one input (bus0)
and/or output (i.e. bus1, bus2, bus3) with a given efficieny (efficiency, efficiency2, efficiency3).
url = "https://tubcloud.tu-berlin.de/s/Rz4iDEQq8DwBd6R/download/electricity-network.nc"
n = pypsa.Network(url)
INFO:pypsa.io:Retrieving network data from https://tubcloud.tu-berlin.de/s/Rz4iDEQq8DwBd6R/download/electricity-network.nc
INFO:pypsa.io:Imported network electricity-network.nc has buses, carriers, generators, global_constraints, loads, storage_units
n
PyPSA Network
Components:
- Bus: 1
- Carrier: 6
- Generator: 4
- GlobalConstraint: 1
- Load: 1
- StorageUnit: 2
Snapshots: 2190
n.buses.index
Index(['electricity'], dtype='object', name='Bus')
n.generators.index
Index(['OCGT', 'onwind', 'offwind', 'solar'], dtype='object', name='Generator')
n.storage_units.index
Index(['battery storage', 'hydrogen storage underground'], dtype='object', name='StorageUnit')
Hydrogen Production#
The following example shows how to model the components of hydrogen storage separately, i.e. electrolysis, fuel cell and storage.
First, let’s remove the simplified hydrogen storage representation:
n.remove("StorageUnit", "hydrogen storage underground")
Add a separate Bus for the hydrogen energy carrier:
n.add("Bus", "hydrogen", carrier='hydrogen')
Index(['hydrogen'], dtype='object')
Add a Link for the hydrogen electrolysis:
n.add(
"Link",
"electrolysis",
bus0="electricity",
bus1="hydrogen",
carrier="electrolysis",
p_nom_extendable=True,
efficiency=0.7,
capital_cost=50e3, # €/MW/a
)
Index(['electrolysis'], dtype='object')
Add a Link for the fuel cell which reconverts hydrogen to electricity:
n.add(
"Link",
"fuel cell",
bus0="hydrogen",
bus1="electricity",
carrier="fuel cell",
p_nom_extendable=True,
efficiency=0.5,
capital_cost=120e3, # €/MW/a
)
Index(['fuel cell'], dtype='object')
Add a Store for the hydrogen storage:
n.add(
"Store",
"hydrogen storage",
bus="hydrogen",
carrier="hydrogen storage",
capital_cost=140, # €/MWh/a
e_nom_extendable=True,
e_cyclic=True, # cyclic state of charge
)
Index(['hydrogen storage'], dtype='object')
We can also add a hydrogen demand to the hydrogen bus.
In the example below, we add a constant hydrogen demand the size of the electricity demand.
p_set = n.loads_t.p_set["demand"].mean()
p_set
54671.88812785388
n.add("Load", "hydrogen demand", bus="hydrogen", carrier="hydrogen", p_set=p_set) # MW
Index(['hydrogen demand'], dtype='object')
When we now optimize the model with additional hydrogen demand…
n.optimize()
WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['electrolysis', 'fuel cell'], dtype='object', name='Link')
WARNING:pypsa.consistency:The following stores have carriers which are not defined:
Index(['hydrogen storage'], dtype='object', name='Store')
WARNING:pypsa.consistency:The following loads have carriers which are not defined:
Index(['hydrogen demand'], dtype='object', name='Load')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity', 'hydrogen'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['electrolysis', 'fuel cell'], dtype='object', name='Link')
WARNING:pypsa.consistency:The following stores have carriers which are not defined:
Index(['hydrogen storage'], dtype='object', name='Store')
WARNING:pypsa.consistency:The following loads have carriers which are not defined:
Index(['hydrogen demand'], dtype='object', name='Load')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity', 'hydrogen'], dtype='object', name='Bus')
/opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/linopy/common.py:147: UserWarning:
coords for dimension(s) ['Generator'] is not aligned with the pandas object. Previously, the indexes of the pandas were ignored and overwritten in these cases. Now, the pandas object's coordinates are taken considered for alignment.
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 0%| | 0/26 [00:00<?, ?it/s]
Writing constraints.: 38%|███▊ | 10/26 [00:00<00:00, 92.71it/s]
Writing constraints.: 77%|███████▋ | 20/26 [00:00<00:00, 76.21it/s]
Writing constraints.: 100%|██████████| 26/26 [00:00<00:00, 66.08it/s]
Writing continuous variables.: 0%| | 0/11 [00:00<?, ?it/s]
Writing continuous variables.: 100%|██████████| 11/11 [00:00<00:00, 160.08it/s]
INFO:linopy.io: Writing time: 0.49s
INFO:linopy.solvers:Log file at /tmp/highs.log
Running HiGHS 1.8.0 (git hash: eda5cbe): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
Matrix [1e-04, 6e+00]
Cost [4e-02, 2e+05]
Bound [0e+00, 0e+00]
RHS [3e+04, 8e+04]
Presolving model
27420 rows, 20857 cols, 73500 nonzeros 0s
25230 rows, 18667 cols, 69120 nonzeros 0s
Presolve : Reductions: rows 25230(-27343); columns 18667(-5431); elements 69120(-45908)
Solving the presolved LP
Using EKK dual simplex solver - serial
Iteration Objective Infeasibilities num(sum)
0 0.0000000000e+00 Pr: 4380(8.00224e+09) 0s
17435 7.0371994395e+10 Pr: 8886(5.10211e+11); Du: 0(1.21526e-08) 5s
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 24098 primals, 52573 duals
Objective: 9.45e+10
Solver model: available
Solver message: optimal
26444 9.4542662147e+10 Pr: 0(0); Du: 0(2.40086e-14) 10s
Solving the original LP from the solution after postsolve
Model status : Optimal
Simplex iterations: 26444
Objective value : 9.4542662147e+10
HiGHS run time : 9.96
Writing the solution to /tmp/linopy-solve-dcjh71u4.sol
INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-ext-p-lower, Generator-ext-p-upper, Link-ext-p-lower, Link-ext-p-upper, Store-ext-e-lower, Store-ext-e-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance, Store-energy_balance were not assigned to the network.
('ok', 'optimal')
… we get an additional price time series for the hydrogen “market”,
n.buses_t.marginal_price.plot()
and can see the individual sizing of the electrolyser, fuel cell and hydrogen storage:
n.statistics.expanded_capacity().div(1e3).round(1)
component carrier
Link electrolysis 183.3
fuel cell 71.3
Store hydrogen storage 41979.1
StorageUnit battery storage 65.1
Generator offwind 179.2
onwind 80.3
solar 480.3
dtype: float64
Furthermore, we might want to explore the storage state of charge of the hydrogen storage and the balancing patterns:
n.stores_t.e.plot()
Heat Demand#
For modelling simple heating systems, we create another bus and connect a load with the heat demand time series to it:
n.add("Bus", "heat", carrier='heat')
Index(['heat'], dtype='object')
url = "https://tubcloud.tu-berlin.de/s/mSkHERH8fJCKNXx/download/heat-load-example.csv"
p_set = pd.read_csv(url, index_col=0, parse_dates=True).squeeze()
p_set.head()
snapshot
2015-01-01 00:00:00 61726.043437
2015-01-01 04:00:00 108787.133591
2015-01-01 08:00:00 101508.988082
2015-01-01 12:00:00 90475.260586
2015-01-01 16:00:00 96307.755312
Name: 0, dtype: float64
n.add("Load", "heat demand", carrier="heat", bus="heat", p_set=p_set)
Index(['heat demand'], dtype='object')
n.loads_t.p_set.div(1e3).plot()
What is now missing are a few heat supply options…
Heat pumps#
To model heat pumps, first we have to calculate the coefficient of performance (COP) profile based on the temperature profile of the heat source.
In the example below, we calculate the COP for an air-sourced heat pump with a sink temperature of 55° C and a population-weighted ambient temperature profile for Germany.
The heat pump performance is assumed to be given by the following function:
where \(\Delta T = T_{sink} - T_{source}\).
def cop(t_source, t_sink=55):
delta_t = t_sink - t_source
return 6.81 - 0.121 * delta_t + 0.000630 * delta_t**2
url = "https://tubcloud.tu-berlin.de/s/S4jRAQMP5Te96jW/download/ninja_weather_country_DE_merra-2_population_weighted.csv"
temp = pd.read_csv(url, skiprows=2, index_col=0, parse_dates=True).loc[
"2015", "temperature"
][::4]
cop(temp).plot()
plt.scatter(temp, cop(temp))
plt.xlabel("temperature [°C]")
plt.ylabel("COP [-]")
Text(0, 0.5, 'COP [-]')
Once we have calculated the heat pump coefficient of performance, we can add the heat pump to the network as a Link. We use the parameter efficiency to incorporate the COP.
n.add(
"Link",
"heat pump",
carrier="heat pump",
bus0="electricity",
bus1="heat",
efficiency=cop(temp),
p_nom_extendable=True,
capital_cost=3e5, # €/MWe/a
)
Index(['heat pump'], dtype='object')
Let’s also add a resistive heater as backup technology:
n.add(
"Link",
"resistive heater",
carrier="resistive heater",
bus0="electricity",
bus1="heat",
efficiency=0.9,
capital_cost=1e4, # €/MWe/a
p_nom_extendable=True,
)
Index(['resistive heater'], dtype='object')
Combined Heat-and-Power (CHP)#
In the following, we are going to add gas-fired combined heat-and-power plants (CHPs). Today, these would use fossil gas, but in the example below we assume imported green methane or methanol with relatively high marginal costs. Since we have no other net emission technology, we can remove the CO\(_2\) limit.
n.remove("GlobalConstraint", "CO2Limit")
Then, we explicitly represent the energy carrier gas:
n.add("Bus", "gas", carrier="gas")
Index(['gas'], dtype='object')
And add a Store of gas, which can be depleted (up to 100 TWh) with fuel costs of 180 €/MWh.
n.add(
"Store",
"gas storage",
carrier="gas storage",
e_initial=100e6, # MWh
e_nom=100e6, # MWh
bus="gas",
marginal_cost=180, # €/MWh_th
)
Index(['gas storage'], dtype='object')
When we do this, we have to model the OCGT power plant as link which converts gas to electricity, not as generator.
n.remove("Generator", "OCGT")
n.add(
"Link",
"OCGT",
bus0="gas",
bus1="electricity",
carrier="OCGT",
p_nom_extendable=True,
capital_cost=20000, # €/MW/a
efficiency=0.4,
)
Index(['OCGT'], dtype='object')
Next, we are going to add a combined heat-and-power (CHP) plant with fixed heat-power ratio (i.e. backpressure operation).
Note
If you want to model flexible heat-power ratios, have a look at this example: https://pypsa.readthedocs.io/en/latest/examples/power-to-gas-boiler-chp.html
n.add(
"Link",
"CHP",
bus0="gas",
bus1="electricity",
bus2="heat",
carrier="CHP",
p_nom_extendable=True,
capital_cost=40000,
efficiency=0.4,
efficiency2=0.4,
)
Index(['CHP'], dtype='object')
Now, let’s optimize the current status of model:
n.optimize()
WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['electrolysis', 'fuel cell', 'heat pump', 'resistive heater', 'CHP'], dtype='object', name='Link')
WARNING:pypsa.consistency:Encountered nan's in static data for columns ['efficiency2'] of component 'Link'.
WARNING:pypsa.consistency:The following stores have carriers which are not defined:
Index(['hydrogen storage', 'gas storage'], dtype='object', name='Store')
WARNING:pypsa.consistency:The following loads have carriers which are not defined:
Index(['hydrogen demand', 'heat demand'], dtype='object', name='Load')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity', 'hydrogen', 'heat', 'gas'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['electrolysis', 'fuel cell', 'heat pump', 'resistive heater', 'CHP'], dtype='object', name='Link')
WARNING:pypsa.consistency:Encountered nan's in static data for columns ['efficiency2'] of component 'Link'.
WARNING:pypsa.consistency:The following stores have carriers which are not defined:
Index(['hydrogen storage', 'gas storage'], dtype='object', name='Store')
WARNING:pypsa.consistency:The following loads have carriers which are not defined:
Index(['hydrogen demand', 'heat demand'], dtype='object', name='Load')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity', 'hydrogen', 'heat', 'gas'], dtype='object', name='Bus')
/opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/linopy/common.py:147: UserWarning:
coords for dimension(s) ['Generator'] is not aligned with the pandas object. Previously, the indexes of the pandas were ignored and overwritten in these cases. Now, the pandas object's coordinates are taken considered for alignment.
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 0%| | 0/27 [00:00<?, ?it/s]
Writing constraints.: 41%|████ | 11/27 [00:00<00:00, 100.28it/s]
Writing constraints.: 81%|████████▏ | 22/27 [00:00<00:00, 79.38it/s]
Writing constraints.: 100%|██████████| 27/27 [00:00<00:00, 59.41it/s]
Writing continuous variables.: 0%| | 0/11 [00:00<?, ?it/s]
Writing continuous variables.: 100%|██████████| 11/11 [00:00<00:00, 150.88it/s]
INFO:linopy.io: Writing time: 0.55s
INFO:linopy.solvers:Log file at /tmp/highs.log
Running HiGHS 1.8.0 (git hash: eda5cbe): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
Matrix [1e-04, 6e+00]
Cost [4e-02, 3e+05]
Bound [0e+00, 0e+00]
RHS [2e+04, 1e+08]
Presolving model
42749 rows, 34000 cols, 119487 nonzeros 0s
38370 rows, 29621 cols, 110729 nonzeros 0s
Presolve : Reductions: rows 38370(-38294); columns 29621(-5430); elements 110729(-50291)
Solving the presolved LP
Using EKK dual simplex solver - serial
Iteration Objective Infeasibilities num(sum)
0 0.0000000000e+00 Ph1: 0(0) 0s
15570 2.2962917512e+10 Pr: 7061(9.5299e+11); Du: 0(1.51356e-06) 5s
21609 3.7705371248e+10 Pr: 8310(1.11402e+13); Du: 0(2.50085e-07) 10s
26715 9.9340691224e+10 Pr: 12410(2.58287e+13); Du: 0(7.68579e-08) 16s
31162 1.3762022038e+11 Pr: 15176(1.90995e+12) 22s
33878 1.4301664789e+11 Pr: 0(0); Du: 0(3.73035e-14) 25s
Solving the original LP from the solution after postsolve
Model status : Optimal
Simplex iterations: 33878
Objective value : 1.4301664789e+11
HiGHS run time : 24.83
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 35051 primals, 76664 duals
Objective: 1.43e+11
Solver model: available
Solver message: optimal
INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-ext-p-lower, Generator-ext-p-upper, Link-ext-p-lower, Link-ext-p-upper, Store-fix-e-lower, Store-fix-e-upper, Store-ext-e-lower, Store-ext-e-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance, Store-energy_balance were not assigned to the network.
Writing the solution to /tmp/linopy-solve-wn7t42yh.sol
('ok', 'optimal')
The objective cost in bn€/a:
n.objective / 1e9
143.01664789325199
The heat energy balance (positive is supply, negative is consumption):
n.statistics.energy_balance(bus_carrier='heat').div(1e6).round(1)
component carrier bus_carrier
Link heat pump heat 519.1
resistive heater heat 34.6
CHP heat 40.0
Load heat heat -593.7
dtype: float64
The heat energy balance as a time series:
n.statistics.energy_balance(aggregate_time=False, bus_carrier='heat').div(1e3).groupby("carrier").sum().T.plot()
The price time series:
n.buses_t.marginal_price.plot()
Long-duration heat storage#
One technology of particular interest in district heating systems with large shares of renewables is long-duration thermal energy storage.
In the following, we are going to introduce a heat storage with investment cost of approximately 3 €/kWh. The energy is not perfectly stored in water tanks. There are standing losses. The decay of thermal energy in the heat storage is modelled through the function \(1 - e^{-\frac{1}{24\tau}}\), where \(\tau\) is assumed to be 180 days. We want to see how that influences the optimal design decisions in the heating sector.
n.add(
"Store",
"heat storage",
bus="heat",
carrier="heat storage",
capital_cost=300, # roughly annuity of 3 €/kWh
standing_loss=1 - np.exp(-1 / 24 / 180),
e_nom_extendable=True,
)
Index(['heat storage'], dtype='object')
n.optimize()
WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['electrolysis', 'fuel cell', 'heat pump', 'resistive heater', 'CHP'], dtype='object', name='Link')
WARNING:pypsa.consistency:Encountered nan's in static data for columns ['efficiency2'] of component 'Link'.
WARNING:pypsa.consistency:The following stores have carriers which are not defined:
Index(['hydrogen storage', 'gas storage', 'heat storage'], dtype='object', name='Store')
WARNING:pypsa.consistency:The following loads have carriers which are not defined:
Index(['hydrogen demand', 'heat demand'], dtype='object', name='Load')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity', 'hydrogen', 'heat', 'gas'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['electrolysis', 'fuel cell', 'heat pump', 'resistive heater', 'CHP'], dtype='object', name='Link')
WARNING:pypsa.consistency:Encountered nan's in static data for columns ['efficiency2'] of component 'Link'.
WARNING:pypsa.consistency:The following stores have carriers which are not defined:
Index(['hydrogen storage', 'gas storage', 'heat storage'], dtype='object', name='Store')
WARNING:pypsa.consistency:The following loads have carriers which are not defined:
Index(['hydrogen demand', 'heat demand'], dtype='object', name='Load')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity', 'hydrogen', 'heat', 'gas'], dtype='object', name='Bus')
/opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/linopy/common.py:147: UserWarning:
coords for dimension(s) ['Generator'] is not aligned with the pandas object. Previously, the indexes of the pandas were ignored and overwritten in these cases. Now, the pandas object's coordinates are taken considered for alignment.
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 0%| | 0/27 [00:00<?, ?it/s]
Writing constraints.: 41%|████ | 11/27 [00:00<00:00, 101.35it/s]
Writing constraints.: 81%|████████▏ | 22/27 [00:00<00:00, 73.57it/s]
Writing constraints.: 100%|██████████| 27/27 [00:00<00:00, 59.19it/s]
Writing continuous variables.: 0%| | 0/11 [00:00<?, ?it/s]
Writing continuous variables.: 100%|██████████| 11/11 [00:00<00:00, 137.22it/s]
INFO:linopy.io: Writing time: 0.56s
INFO:linopy.solvers:Log file at /tmp/highs.log
Running HiGHS 1.8.0 (git hash: eda5cbe): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
Matrix [1e-04, 6e+00]
Cost [4e-02, 3e+05]
Bound [0e+00, 0e+00]
RHS [2e+04, 1e+08]
Presolving model
47128 rows, 38380 cols, 132624 nonzeros 0s
40560 rows, 31812 cols, 119488 nonzeros 0s
40560 rows, 31812 cols, 119488 nonzeros 0s
Presolve : Reductions: rows 40560(-42675); columns 31812(-7620); elements 119488(-56862)
Solving the presolved LP
Using EKK dual simplex solver - serial
Iteration Objective Infeasibilities num(sum)
0 0.0000000000e+00 Ph1: 0(0) 0s
16443 7.6508187661e+09 Pr: 8788(1.48947e+14); Du: 0(9.70147e-07) 5s
22711 4.3725291538e+10 Pr: 6314(1.40186e+13); Du: 0(1.2653e-06) 11s
27869 4.4056163819e+10 Pr: 7072(4.98632e+14); Du: 0(6.05973e-07) 16s
31942 6.7369736899e+10 Pr: 10332(3.3962e+15); Du: 0(2.37077e-07) 21s
35534 9.5675850589e+10 Pr: 18650(1.00973e+14) 27s
The objective cost (in bn€/a) was reduced by 7%!
n.objective / 1e9
133.41557805470643
The heat energy balance shows the additional losses of the heat storage and the added supply:
n.statistics.energy_balance(bus_carrier='heat').div(1e6).round(1)
component carrier bus_carrier
Load heat heat -593.7
Store heat storage heat -11.7
Link heat pump heat 499.2
resistive heater heat 82.3
CHP heat 23.8
dtype: float64
The heat energy balance as a time series:
n.statistics.energy_balance(aggregate_time=False, bus_carrier='heat').div(1e3).groupby("carrier").sum().T.plot()
The heat price time series is much less impacted by the variability of the electricity price:
n.buses_t.marginal_price.plot()
The different storage state of charge time series:
n.stores_t.e.plot()
Electric Vehicles#
To model electric vehicles, we first create another bus for the electric vehicles.
n.add("Bus", "EV", carrier="EV")
Index(['EV'], dtype='object')
Then, we can attach the electricity consumption of electric vehicles to this bus:
url = "https://tubcloud.tu-berlin.de/s/9r5bMSbzzQiqG7H/download/electric-vehicle-profile-example.csv"
p_set = pd.read_csv(url, index_col=0, parse_dates=True).squeeze()
p_set.loc["2015-01"].div(1e3).plot()
n.add("Load", "EV demand", bus="EV", carrier="EV demand", p_set=p_set)
Index(['EV demand'], dtype='object')
The electric vehicles can only be charged when they are plugged-in. Below we load an availability profile telling us what share of electric vehicles is plugged-in at home – we only assume home charging in this example.
url = "https://tubcloud.tu-berlin.de/s/E3PBWPfYaWwCq7a/download/electric-vehicle-availability-example.csv"
availability_profile = pd.read_csv(url, index_col=0, parse_dates=True).squeeze()
availability_profile.loc["2015-01"].plot()
Then, we can add a link for the electric vehicle charger using assumption about the number of EVs and their charging rates.
number_cars = 40e6 # number of EV cars
bev_charger_rate = 0.011 # 3-phase EV charger with 11 kW
p_nom = number_cars * bev_charger_rate
n.add(
"Link",
"EV charger",
bus0="electricity",
bus1="EV",
p_nom=p_nom,
carrier="EV charger",
p_max_pu=availability_profile,
efficiency=0.9,
)
Index(['EV charger'], dtype='object')
We can also allow vehicle-to-grid operation (i.e. electric vehicles inject power into the grid):
n.add(
"Link",
"V2G",
bus0="EV",
bus1="electricity",
p_nom=p_nom,
carrier="V2G",
p_max_pu=availability_profile,
efficiency=0.9,
)
Index(['V2G'], dtype='object')
The demand-side management potential we model as a store. This is not unlike a battery storage, but we impose additional constraints on when the store needs to be charged to a certain level (e.g. 75% full every morning).
bev_energy = 0.05 # average battery size of EV in MWh
bev_dsm_participants = 0.5 # share of cars that do smart charging
e_nom = number_cars * bev_energy * bev_dsm_participants
url = "https://tubcloud.tu-berlin.de/s/K62yACBRTrxLTia/download/dsm-profile-example.csv"
dsm_profile = pd.read_csv(url, index_col=0, parse_dates=True).squeeze().shift(2).fillna(0)[::4]
dsm_profile.loc["2015-01-01":"2015-01-07"].plot()
n.add(
"Store",
"EV DSM",
bus="EV",
carrier="EV battery",
e_cyclic=True, # state of charge at beginning = state of charge at the end
e_nom=e_nom,
e_min_pu=dsm_profile,
)
Index(['EV DSM'], dtype='object')
Then, we can solve the fully sector-coupled model altogether including electricity, passenger transport, hydrogen and heating.
n.optimize(solver_name="highs")
Show code cell output
WARNING:pypsa.consistency:The following stores have carriers which are not defined:
Index(['hydrogen storage', 'gas storage', 'heat storage', 'EV DSM'], dtype='object', name='Store')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity', 'hydrogen', 'heat', 'gas', 'EV'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following loads have carriers which are not defined:
Index(['hydrogen demand', 'heat demand', 'EV demand'], dtype='object', name='Load')
WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['electrolysis', 'fuel cell', 'heat pump', 'resistive heater', 'CHP',
'EV charger', 'V2G'],
dtype='object', name='Link')
WARNING:pypsa.consistency:The following stores have carriers which are not defined:
Index(['hydrogen storage', 'gas storage', 'heat storage', 'EV DSM'], dtype='object', name='Store')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity', 'hydrogen', 'heat', 'gas', 'EV'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following loads have carriers which are not defined:
Index(['hydrogen demand', 'heat demand', 'EV demand'], dtype='object', name='Load')
WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['electrolysis', 'fuel cell', 'heat pump', 'resistive heater', 'CHP',
'EV charger', 'V2G'],
dtype='object', name='Link')
/home/fneum/miniconda3/envs/nordnet/lib/python3.12/site-packages/linopy/common.py:147: UserWarning:
coords for dimension(s) ['Generator'] is not aligned with the pandas object. Previously, the indexes of the pandas were ignored and overwritten in these cases. Now, the pandas object's coordinates are taken considered for alignment.
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 29/29 [00:00<00:00, 84.81it/s]
Writing continuous variables.: 100%|██████████| 11/11 [00:00<00:00, 226.05it/s]
INFO:linopy.io: Writing time: 0.41s
INFO:linopy.solvers:Log file at /tmp/highs.log
Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
Matrix [1e-04, 6e+00]
Cost [4e-02, 3e+05]
Bound [0e+00, 0e+00]
RHS [2e+03, 1e+08]
Presolving model
51508 rows, 47140 cols, 150144 nonzeros 0s
42750 rows, 38382 cols, 132628 nonzeros 0s
42750 rows, 38382 cols, 132628 nonzeros 0s
Presolve : Reductions: rows 42750(-58005); columns 38382(-9810); elements 132628(-74382)
Solving the presolved LP
Using EKK dual simplex solver - serial
Iteration Objective Infeasibilities num(sum)
0 0.0000000000e+00 Ph1: 0(0) 0s
20489 7.0209292127e+09 Pr: 10189(7.06301e+14); Du: 0(9.85806e-07) 5s
26917 3.1666331621e+10 Pr: 6630(2.35325e+13); Du: 0(6.23236e-07) 11s
32263 3.2859044901e+10 Pr: 7443(6.50205e+13); Du: 0(2.8441e-07) 16s
36620 5.3585544208e+10 Pr: 10737(1.8223e+13); Du: 0(4.035e-07) 21s
40843 8.8898698525e+10 Pr: 15887(1.50189e+15); Du: 0(7.42276e-07) 26s
45777 1.0089919205e+11 Pr: 9171(1.94314e+16); Du: 0(5.04967e-07) 32s
49643 1.1413934646e+11 Pr: 14967(2.15477e+13); Du: 0(5.47685e-07) 37s
53329 1.2988564883e+11 Pr: 10975(1.23989e+13); Du: 0(6.79277e-07) 43s
56786 1.3519466178e+11 Pr: 4851(3.46768e+11); Du: 0(5.56497e-07) 48s
59796 1.3651637401e+11 Pr: 0(0); Du: 0(2.27374e-13) 52s
59796 1.3651637401e+11 Pr: 0(0); Du: 0(2.27374e-13) 52s
Solving the original LP from the solution after postsolve
Model status : Optimal
Simplex iterations: 59796
Objective value : 1.3651637401e+11
HiGHS run time : 52.41
Writing the solution to /tmp/linopy-solve-tk463w7g.sol
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 48192 primals, 100755 duals
Objective: 1.37e+11
Solver model: available
Solver message: optimal
INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-ext-p-lower, Generator-ext-p-upper, Link-fix-p-lower, Link-fix-p-upper, Link-ext-p-lower, Link-ext-p-upper, Store-fix-e-lower, Store-fix-e-upper, Store-ext-e-lower, Store-ext-e-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance, Store-energy_balance were not assigned to the network.
('ok', 'optimal')
A very general overview is given with a call to n.statistics():
n.statistics()
| Optimal Capacity | Installed Capacity | Supply | Withdrawal | Energy Balance | Transmission | Capacity Factor | Curtailment | Capital Expenditure | Operational Expenditure | Revenue | Market Value | ||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Generator | offwind | 2.455913e+05 | 0.0 | 7.805015e+08 | 0.000000e+00 | 7.805015e+08 | 0.0 | 0.362791 | 0.000000e+00 | 4.286949e+10 | 1.654663e+07 | 4.288603e+10 | 54.946768 |
| onwind | 1.575662e+05 | 0.0 | 2.659010e+08 | 0.000000e+00 | 2.659010e+08 | 0.0 | 0.192643 | 1.808124e+07 | 1.601567e+10 | 3.798662e+08 | 1.639554e+10 | 61.660318 | |
| solar | 5.033996e+05 | 0.0 | 5.476326e+08 | 0.000000e+00 | 5.476326e+08 | 0.0 | 0.124186 | 0.000000e+00 | 2.584798e+10 | 5.804905e+06 | 2.585378e+10 | 47.210086 | |
| Link | CHP | 1.123002e+05 | 0.0 | 5.083544e+07 | 6.354431e+07 | -1.270886e+07 | 0.0 | 0.064594 | 0.000000e+00 | 4.492008e+09 | 0.000000e+00 | 4.492008e+09 | 88.363693 |
| EV charger | 4.400000e+05 | 440000.0 | 2.000217e+08 | 2.222463e+08 | -2.222463e+07 | 0.0 | 0.057660 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | NaN | |
| OCGT | 6.770876e+03 | 0.0 | 3.536607e+05 | 8.841518e+05 | -5.304911e+05 | 0.0 | 0.014907 | 0.000000e+00 | 1.354175e+08 | 0.000000e+00 | 1.354175e+08 | 382.902358 | |
| V2G | 4.400000e+05 | 440000.0 | 6.816967e+07 | 7.574408e+07 | -7.574408e+06 | 0.0 | 0.019651 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | NaN | |
| electrolysis | 1.579354e+05 | 0.0 | 4.789257e+08 | 6.841796e+08 | -2.052539e+08 | 0.0 | 0.494523 | 0.000000e+00 | 7.896770e+09 | 0.000000e+00 | 7.896770e+09 | 16.488505 | |
| heat pump | 4.444094e+04 | 0.0 | 4.804672e+08 | 1.924107e+08 | 2.880565e+08 | 0.0 | 0.494244 | 0.000000e+00 | 1.333228e+10 | 0.000000e+00 | 1.333228e+10 | 27.748578 | |
| resistive heater | 9.724606e+04 | 0.0 | 9.919236e+07 | 1.102137e+08 | -1.102137e+07 | 0.0 | 0.129378 | 0.000000e+00 | 9.724606e+08 | 0.000000e+00 | 9.724606e+08 | 9.803785 | |
| Load | - | 0.000000e+00 | 0.0 | 0.000000e+00 | 4.789257e+08 | -4.789257e+08 | 0.0 | NaN | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | -4.218900e+10 | NaN |
| EV demand | 0.000000e+00 | 0.0 | 0.000000e+00 | 1.242776e+08 | -1.242776e+08 | 0.0 | NaN | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | -1.104624e+10 | NaN | |
| heat | 0.000000e+00 | 0.0 | 0.000000e+00 | 5.936692e+08 | -5.936692e+08 | 0.0 | NaN | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | -4.801251e+10 | NaN | |
| hydrogen | 0.000000e+00 | 0.0 | 0.000000e+00 | 4.789257e+08 | -4.789257e+08 | 0.0 | NaN | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | -3.838917e+10 | NaN | |
| Store | EV battery | 1.000000e+06 | 1000000.0 | 1.290614e+08 | 1.290614e+08 | 0.000000e+00 | 0.0 | 0.726787 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 3.120540e+09 | 24.178714 |
| gas storage | 1.000000e+08 | 100000000.0 | 6.442846e+07 | 0.000000e+00 | 6.442846e+07 | 0.0 | 0.499883 | 0.000000e+00 | 0.000000e+00 | 1.159712e+10 | 1.159712e+10 | 180.000000 | |
| heat storage | 1.270014e+07 | 0.0 | 1.953699e+08 | 2.067779e+08 | -1.140802e+07 | 0.0 | 0.443182 | 0.000000e+00 | 3.810042e+09 | 0.000000e+00 | 3.810042e+09 | 19.501679 | |
| hydrogen storage | 6.532086e+07 | 0.0 | 2.009489e+08 | 2.009489e+08 | 0.000000e+00 | 0.0 | 0.565556 | 0.000000e+00 | 9.144920e+09 | 0.000000e+00 | 9.144920e+09 | 45.508680 |
And it can tell you statistics about the capital expenditures:
n.statistics.capex().groupby("carrier").sum().div(1e9).sort_values().dropna().plot.bar()
Exercises#
Explore how the model reacts to changing assumptions and available technologies. Here are a few inspirations, but choose in any order according to your interests:
By how much did the optimised capacities react to
Assume underground hydrogen storage is not geographically available. Increase the cost of hydrogen storage by factor 10. How does the model react?
Add a ground-sourced heat pump with a constant COP function of 3.5 but double the investment costs. Would this technology get built? How low would the costs need to be?
Vary the charging rates and availability profiles of electric vehicles, and whether they participate in demand-side management (V2G, smart charging) or not. How does it affect the system design?
Limit green gas imports to 10 TWh. What does the model do in periods with persistent low wind and solar feed-in but high heating demand?
You can also combine these aspects with investigations suggested for the electricity-only model.